Take-Home Exercise 4

Author

Imran Ibrahim

Published

February 25, 2024

Modified

March 1, 2024

Getting Started

Loading R package

pacman::p_load(sf, tidyverse, tmap, dplyr,
               raster, spatstat, spdep,
               lubridate, leaflet,
               plotly, DT, viridis,
               ggplot2, sfdep)

Importing the ACLED data

Country specific data from the Armed Conflict Location & Event Data Project (ACLED) can be downloaded at https://acleddata.com/data-export-tool/

ACLED_MMR <- read_csv("data/MMR.csv")
class(ACLED_MMR)
[1] "spec_tbl_df" "tbl_df"      "tbl"         "data.frame" 

Downloading and loading the shape files for country

Shape files can be downloaded from the Myanmmar Information Management Unit (MIMU) website at https://geonode.themimu.info/layers/?limit=100&offset=0

This source was chosen over GADM and GeoBoundaries due to its updated administrative region information and map levels.

Important- Data Quality Issue with ACLED data

ACLED captures event data from national, sub-national and other credible media sources, and populates event locations based on the last known information.

However, due to the dynamic nature of conflict and politics, country/administrative boundaries and borders can sometimes be fluid. Names of administrative areas were found to have changed; either disaggregated into new countries/administrative areas or previously active but now defunct. Further, some administrative areas were agglomerated and upgraded into higher tier administrative areas.

As part of our data cleaning and preparation process, I had to identify discrepancies in both admin1 and admin2 data levels and re-name some administrative areas to sync with my downloaded shape files from MIMU.

Data Preparation and Cleaning

Loading Admin1(administrative region/area) shape files

mmr_shp_mimu_1 <-  st_read(dsn = "data/geospatial3",  
                  layer = "mmr_polbnda2_adm1_250k_mimu_1")
Reading layer `mmr_polbnda2_adm1_250k_mimu_1' from data source 
  `C:\imranmi\ISSS608-VAA\Take-home-ex\Take-home-Ex4\data\geospatial3' 
  using driver `ESRI Shapefile'
Simple feature collection with 18 features and 6 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 92.1721 ymin: 9.696844 xmax: 101.17 ymax: 28.54554
Geodetic CRS:  WGS 84
class(mmr_shp_mimu_1)
[1] "sf"         "data.frame"

The Shape file for admin1 level map, is an SF object, with geometry type: Multipolygon

st_geometry(mmr_shp_mimu_1)
Geometry set for 18 features 
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 92.1721 ymin: 9.696844 xmax: 101.17 ymax: 28.54554
Geodetic CRS:  WGS 84
First 5 geometries:
unique_regions_mimu1 <- unique(mmr_shp_mimu_1$ST)

unique_regions_mimu1
 [1] "Ayeyarwady"   "Bago (East)"  "Bago (West)"  "Chin"         "Kachin"      
 [6] "Kayah"        "Kayin"        "Magway"       "Mandalay"     "Mon"         
[11] "Nay Pyi Taw"  "Rakhine"      "Sagaing"      "Shan (East)"  "Shan (North)"
[16] "Shan (South)" "Tanintharyi"  "Yangon"      

There are 18 admin1 levels or states/regions in mmr_shp_mimu_1

Lets compare with our admin1 levels in our main dataset ACLED_MMR

unique_acled_regions1 <- unique(ACLED_MMR$admin1)

unique_acled_regions1
 [1] "Rakhine"     "Bago-East"   "Sagaing"     "Shan-North"  "Mandalay"   
 [6] "Mon"         "Yangon"      "Shan-South"  "Kayin"       "Kachin"     
[11] "Magway"      "Ayeyarwady"  "Nay Pyi Taw" "Kayah"       "Chin"       
[16] "Bago-West"   "Tanintharyi" "Shan-East"  

I will write a simple function to identify the discrepancies between the shape file and the region names in our main dataset.

# Find the unique region names that are in 'unique_acled_regions1' but not in 'unique_regions_mimu1'

mismatched_admin1 <- setdiff(unique_acled_regions1, unique_regions_mimu1)

if (length(mismatched_admin1) > 0) {
  print("The following region names from 'acled_mmr' do not match any in 'mimu1':")
  print(mismatched_admin1)
} else {
  print("All unique region names in 'acled_mmr' match the unique region names in 'mimu1.'")
}
[1] "The following region names from 'acled_mmr' do not match any in 'mimu1':"
[1] "Bago-East"  "Shan-North" "Shan-South" "Bago-West"  "Shan-East" 

Lets harmonize the names in both data files. I will resave it to a new data set called ACLED_MMR_1

ACLED_MMR_1 <- ACLED_MMR %>%
  mutate(admin1 = case_when(
    admin1 == "Bago-East" ~ "Bago (East)",
    admin1 == "Bago-West" ~ "Bago (West)",
    admin1 == "Shan-North" ~ "Shan (North)",
    admin1 == "Shan-South" ~ "Shan (South)",
    admin1 == "Shan-East" ~ "Shan (East)",
    TRUE ~ as.character(admin1)
  ))

Checking if our changes are successful.

# Get unique admin 1 region names from 'ACLED_MMR_1'
unique_acled_regions1 <- unique(ACLED_MMR_1$admin1)

# Get unique region names from 'mmr_shp_mimu_1'
unique_map_regions_mimu1 <- unique(mmr_shp_mimu_1$ST)

# Find the unique region names that are in 'unique_acled_regions1' but not in 'unique_map_regions_mimu1'

mismatched_regions <- setdiff(unique_acled_regions1, unique_map_regions_mimu1)

if (length(mismatched_regions) > 0) {
  print("The following region names from 'acled_mmr_1' do not match any in 'mmr_shp_mimu_1':")
  print(mismatched_regions)
} else {
  print("All unique region names in 'acled_mmr_1' match the unique region names in 'mmmr_shp_mimu_1.'")
}
[1] "All unique region names in 'acled_mmr_1' match the unique region names in 'mmmr_shp_mimu_1.'"

Lets do a sample plot to see how our country map looks like at admin1 level

plot(mmr_shp_mimu_1)

Loading Admin2 (administrative region/area) shape files

mmr_shp_mimu_2 <-  st_read(dsn = "data/geospatial3",  
                  layer = "mmr_polbnda_adm2_250k_mimu")
Reading layer `mmr_polbnda_adm2_250k_mimu' from data source 
  `C:\imranmi\ISSS608-VAA\Take-home-ex\Take-home-Ex4\data\geospatial3' 
  using driver `ESRI Shapefile'
Simple feature collection with 80 features and 7 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 92.1721 ymin: 9.696844 xmax: 101.17 ymax: 28.54554
Geodetic CRS:  WGS 84
class(mmr_shp_mimu_2)
[1] "sf"         "data.frame"

The Shape file for admin2 level map, is an SF object, with geometry type: Multipolygon

st_geometry(mmr_shp_mimu_2)
Geometry set for 80 features 
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 92.1721 ymin: 9.696844 xmax: 101.17 ymax: 28.54554
Geodetic CRS:  WGS 84
First 5 geometries:
unique_regions_mimu2 <- unique(mmr_shp_mimu_2$DT)

unique_regions_mimu2
 [1] "Hinthada"                        "Labutta"                        
 [3] "Maubin"                          "Myaungmya"                      
 [5] "Pathein"                         "Pyapon"                         
 [7] "Bago"                            "Taungoo"                        
 [9] "Pyay"                            "Thayarwady"                     
[11] "Falam"                           "Hakha"                          
[13] "Matupi"                          "Mindat"                         
[15] "Bhamo"                           "Mohnyin"                        
[17] "Myitkyina"                       "Puta-O"                         
[19] "Bawlake"                         "Loikaw"                         
[21] "Hpa-An"                          "Hpapun"                         
[23] "Kawkareik"                       "Myawaddy"                       
[25] "Gangaw"                          "Magway"                         
[27] "Minbu"                           "Pakokku"                        
[29] "Thayet"                          "Kyaukse"                        
[31] "Maungdaw"                        "Mrauk-U"                        
[33] "Sittwe"                          "Thandwe"                        
[35] "Hkamti"                          "Kale"                           
[37] "Kanbalu"                         "Katha"                          
[39] "Kawlin"                          "Mawlaik"                        
[41] "Monywa"                          "Naga Self-Administered Zone"    
[43] "Sagaing"                         "Shwebo"                         
[45] "Tamu"                            "Yinmarbin"                      
[47] "Kengtung"                        "Monghsat"                       
[49] "Tachileik"                       "Hopang"                         
[51] "Kokang Self-Administered Zone"   "Kyaukme"                        
[53] "Lashio"                          "Matman"                         
[55] "Mongmit"                         "Muse"                           
[57] "Pa Laung Self-Administered Zone" "Danu Self-Administered Zone"    
[59] "Langkho"                         "Loilen"                         
[61] "Pa-O Self-Administered Zone"     "Taunggyi"                       
[63] "Dawei"                           "Kawthoung"                      
[65] "Mandalay"                        "Meiktila"                       
[67] "Myingyan"                        "Nyaung-U"                       
[69] "Pyinoolwin"                      "Yamethin"                       
[71] "Mawlamyine"                      "Thaton"                         
[73] "Det Khi Na"                      "Oke Ta Ra"                      
[75] "Kyaukpyu"                        "Myeik"                          
[77] "Yangon (East)"                   "Yangon (North)"                 
[79] "Yangon (South)"                  "Yangon (West)"                  

There are 80 admin2 levels or states/districts in mmr_shp_mimu_2

Lets compare with our admin2 levels in our main dataset ACLED_MMR

unique_acled_regions2 <- unique(ACLED_MMR$admin2)

unique_acled_regions2
 [1] "Maungdaw"                        "Bago"                           
 [3] "Shwebo"                          "Kyaukme"                        
 [5] "Pyinoolwin"                      "Muse"                           
 [7] "Sittwe"                          "Yinmarbin"                      
 [9] "Thaton"                          "Yangon-North"                   
[11] "Pa-O Self-Administered Zone"     "Hpapun"                         
[13] "Kyaukpyu"                        "Yangon-West"                    
[15] "Mongmit"                         "Bhamo"                          
[17] "Mrauk-U"                         "Yangon-East"                    
[19] "Yangon-South"                    "Monywa"                         
[21] "Gangaw"                          "Pathein"                        
[23] "Katha"                           "Taungoo"                        
[25] "Kanbalu"                         "Lashio"                         
[27] "Mawlamyine"                      "Myitkyina"                      
[29] "Kawkareik"                       "Loilen"                         
[31] "Mandalay"                        "Kawlin"                         
[33] "Kyaukse"                         "Magway"                         
[35] "Meiktila"                        "Pakokku"                        
[37] "Taunggyi"                        "Tamu"                           
[39] "Nay Pyi Taw"                     "Mohnyin"                        
[41] "Kale"                            "Det Khi Na"                     
[43] "Myingyan"                        "Loikaw"                         
[45] "Matupi"                          "Pyay"                           
[47] "Sagaing"                         "Myeik"                          
[49] "Dawei"                           "Thayarwady"                     
[51] "Thandwe"                         "Mawlaik"                        
[53] "Bawlake"                         "Pyapon"                         
[55] "Hinthada"                        "Thayet"                         
[57] "Pa Laung Self-Administered Zone" "Mindat"                         
[59] "Hkamti"                          "Kokang Self-Administered Zone"  
[61] "Hpa-An"                          "Danu Self-Administered Zone"    
[63] "Myawaddy"                        "Maubin"                         
[65] "Hakha"                           "Falam"                          
[67] "Minbu"                           "Monghsat"                       
[69] "Puta-O"                          "Hopang"                         
[71] "Nyaung-U"                        "Kawthoung"                      
[73] "Yamethin"                        "Yangon"                         
[75] "Myaungmya"                       "Mong Pawk (Wa SAD)"             
[77] "Oke Ta Ra"                       "Matman"                         
[79] "Kengtung"                        "Naga Self-Administered Zone"    
[81] "Labutta"                         "Langkho"                        
[83] "Tachileik"                      

I will write a simple function to identify the discrepancies between the shape file and our state/district names in our main dataset.

# Find the unique region names that are in 'unique_acled_regions2' but not in 'unique_regions_mimu2'

mismatched_admin2 <- setdiff(unique_acled_regions2, unique_regions_mimu2)

if (length(mismatched_admin2) > 0) {
  print("The following region names from 'acled_mmr' do not match any in 'mimu2':")
  print(mismatched_admin2)
} else {
  print("All unique region names in 'acled_mmr' match the unique region names in 'mimu2.'")
}
[1] "The following region names from 'acled_mmr' do not match any in 'mimu2':"
[1] "Yangon-North"       "Yangon-West"        "Yangon-East"       
[4] "Yangon-South"       "Nay Pyi Taw"        "Yangon"            
[7] "Mong Pawk (Wa SAD)"

Lets harmonize the names in both data files. I will resave it to the previous data set called ACLED_MMR_1

ACLED_MMR_1 <- ACLED_MMR_1 %>%
  mutate(admin2 = case_when(
    admin2 == "Yangon-East" ~ "Yangon (East)",
    admin2 == "Yangon-West" ~ "Yangon (West)",
    admin2 == "Yangon-North" ~ "Yangon (North)",
    admin2 == "Yangon-South" ~ "Yangon (South)",
    admin2 == "Mong Pawk (Wa SAD)" ~ "Tachileik",
    admin2 == "Nay Pyi Taw" ~ "Det Khi Na",
    admin2 == "Yangon" ~ "Yangon (West)",
    TRUE ~ as.character(admin2)
  ))

Checking if our changes are successful.

# Get unique admin 2 district names from 'ACLED_MMR_1'
unique_acled_regions2 <- unique(ACLED_MMR_1$admin2)

# Get unique district names from 'mmr_shp_mimu_2'
unique_map_regions_mimu2 <- unique(mmr_shp_mimu_2$DT)

# Find the unique district names that are in 'unique_acled_regions2' but not in 'unique_map_regions_mimu2'

mismatched_regions2 <- setdiff(unique_acled_regions2, unique_map_regions_mimu2)

if (length(mismatched_regions2) > 0) {
  print("The following district names from 'acled_mmr_1' do not match any in 'mmr_shp_mimu_2':")
  print(mismatched_regions2)
} else {
  print("All unique district names in 'acled_mmr_1' match the unique district names in 'mmmr_shp_mimu_2.'")
}
[1] "All unique district names in 'acled_mmr_1' match the unique district names in 'mmmr_shp_mimu_2.'"

Lets do a sample plot to see how our country map looks like at admin2 (districts) level.

plot(mmr_shp_mimu_2)

Data Wrangling

For the purposes of plotting choropleth maps, I will first create subsets for each admin level

data_prepared_1 <- ACLED_MMR_1 %>%
    group_by(year, admin1) %>%
    summarise(Incidents = n(),
              Fatalities = sum(fatalities, na.rm = TRUE)) %>%
              
    ungroup()

Checking the total sums of incidents and fatalities.

# Calculate the total number of incidents
total_incidents1 <- sum(data_prepared_1$Incidents)

# Calculate the total number of fatalities
total_fatalities1 <- sum(data_prepared_1$Fatalities)

# Print the totals
print(paste("Total Incidents:", total_incidents1))
[1] "Total Incidents: 57198"
print(paste("Total Fatalities:", total_fatalities1))
[1] "Total Fatalities: 57593"
data_prepared_2 <- ACLED_MMR_1 %>%
    group_by(year, admin2) %>%
    summarise(Incidents = n(),
              Fatalities = sum(fatalities, na.rm = TRUE)) %>%
              
    ungroup()

Checking the total sums of incidents and fatalities.

# Calculate the total number of incidents
total_incidents2 <- sum(data_prepared_2$Incidents)

# Calculate the total number of fatalities
total_fatalities2 <- sum(data_prepared_2$Fatalities)

# Print the totals
print(paste("Total Incidents:", total_incidents2))
[1] "Total Incidents: 57198"
print(paste("Total Fatalities:", total_fatalities2))
[1] "Total Fatalities: 57593"

Next, I will do a spatial join between my shape files and attribute files

ACLED_MMR_admin1 <- left_join(mmr_shp_mimu_1, data_prepared_1,
                            by = c("ST" = "admin1"))

ACLED_MMR_admin2 <- left_join(mmr_shp_mimu_2, data_prepared_2,
                            by = c("DT" = "admin2"))
class(ACLED_MMR_admin1)
[1] "sf"         "data.frame"
class(ACLED_MMR_admin2)
[1] "sf"         "data.frame"

Checking that the total number of Incidents and Fatalities in the new datasets match the original attribute file

total_incidents3 <- sum(ACLED_MMR_admin1$Incidents)

# Calculate the total number of fatalities
total_fatalities3 <- sum(ACLED_MMR_admin1$Fatalities)

# Print the totals
print(paste("Total Incidents:", total_incidents3))
[1] "Total Incidents: 57198"
print(paste("Total Fatalities:", total_fatalities3))
[1] "Total Fatalities: 57593"
total_incidents4 <- sum(ACLED_MMR_admin2$Incidents)

# Calculate the total number of fatalities
total_fatalities4 <- sum(ACLED_MMR_admin2$Fatalities)

# Print the totals
print(paste("Total Incidents:", total_incidents4))
[1] "Total Incidents: 57198"
print(paste("Total Fatalities:", total_fatalities4))
[1] "Total Fatalities: 57593"

Choropleth Maps

Choropleth Map of Incidents & Fatalities by Admin1 (Region/State)

In the Shiny App, the below choropleth maps can be plotted, with users being able to choose the following:-

  • variable ie Incidents or Fatalities

  • specific year,

  • year range,

  • data classification type, and

  • number of clusters (n)

tm_shape(ACLED_MMR_admin1)+
  tm_fill("Incidents",
          n = 5,
          style = "quantile") +
  tm_borders(alpha = 0.5)

Choropleth map of Incidents & Fatalities by Admin2 level (by District)

Similarly, in the Shiny App, the below choropleth maps can be plotted, with users being able to choose the following:-

  • Variable ie Incidents or Fatalities

  • specific year,

  • year range,

  • data classification type, and

  • number of clusters (n)

tm_shape(ACLED_MMR_admin2)+
  tm_fill("Incidents",
          n = 6,
          style = "quantile",
          palette = "Reds") +
  tm_borders(alpha = 0.5)

tm_shape(ACLED_MMR_admin2)+
  tm_fill("Fatalities",
          n = 6,
          style = "quantile",
          palette = "Reds") +
  tm_borders(alpha = 0.5)

Incidents and Fatalities by individual regions using tm_facet()

tm_shape(ACLED_MMR_admin1) +
  tm_fill("Fatalities",
          style = "quantile",
          palette = "Reds",
          thres.poly = 0) + 
  tm_facets(by="ST", 
            free.coords=TRUE, 
            drop.shapes=TRUE) +
  tm_layout(legend.show = FALSE,
            title.position = c("center", "center"), 
            title.size = 20) +
  tm_borders(alpha = 0.5)

tm_shape(ACLED_MMR_admin1) +
  tm_fill("Incidents",
          style = "quantile",
          palette = "Reds",
          thres.poly = 0) + 
  tm_facets(by="ST", 
            free.coords=TRUE, 
            drop.shapes=TRUE) +
  tm_layout(legend.show = FALSE,
            title.position = c("center", "center"), 
            title.size = 20) +
  tm_borders(alpha = 0.5)

Instead on the plots above, which does not provide specific information on the trends, we can also consider the below line charts.

ggplot(ACLED_MMR_admin1, aes(x = year, y = Fatalities, group = ST, color = ST)) +
  geom_line() + # Creates line plots
  facet_wrap(~ ST, scales = "free_y") + # Creates a separate plot for each region
  theme_minimal() + # Uses a minimal theme for the plot
  labs(title = "Fatalities Over Time by Region",
       x = "Time",
       y = "Fatalities") +
  theme(legend.position = "none")

ggplot(ACLED_MMR_admin1, aes(x = year, y = Incidents, group = ST, color = ST)) +
  geom_line() + # Creates line plots
  facet_wrap(~ ST, scales = "free_y") + # Creates a separate plot for each region
  theme_minimal() + # Uses a minimal theme for the plot
  labs(title = "Incidents Over Time by Region",
       x = "Time",
       y = "Fatalities") +
  theme(legend.position = "none")

Proportional Symbol Maps

Proportional symbol maps (also known as graduate symbol maps) are a class of maps that use the visual variable of size to represent differences in the magnitude of a discrete, abruptly changing phenomenon, e.g. counts of people.

First I will convert the ACLED_MMR_1 data set to become an sf object

# Convert conflict data to an sf object
conflict_sf <- st_as_sf(ACLED_MMR_1, coords = c("longitude", "latitude"), crs = 4326)
class(conflict_sf)
[1] "sf"         "tbl_df"     "tbl"        "data.frame"
conflict_sf
Simple feature collection with 57198 features and 33 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 92.199 ymin: 9.9824 xmax: 100.3576 ymax: 27.731
Geodetic CRS:  WGS 84
# A tibble: 57,198 × 34
   event_id_cnty event_date        year time_precision disorder_type  event_type
 * <chr>         <chr>            <dbl>          <dbl> <chr>          <chr>     
 1 MMR58558      16 February 2024  2024              1 Political vio… Battles   
 2 MMR58559      16 February 2024  2024              1 Political vio… Battles   
 3 MMR58443      16 February 2024  2024              2 Political vio… Violence …
 4 MMR58502      16 February 2024  2024              1 Political vio… Violence …
 5 MMR58507      16 February 2024  2024              1 Political vio… Explosion…
 6 MMR58508      16 February 2024  2024              1 Political vio… Explosion…
 7 MMR58547      16 February 2024  2024              1 Strategic dev… Strategic…
 8 MMR58560      16 February 2024  2024              1 Political vio… Battles   
 9 MMR58589      16 February 2024  2024              2 Political vio… Violence …
10 MMR58648      16 February 2024  2024              1 Political vio… Explosion…
# ℹ 57,188 more rows
# ℹ 28 more variables: sub_event_type <chr>, actor1 <chr>, assoc_actor_1 <chr>,
#   inter1 <dbl>, actor2 <chr>, assoc_actor_2 <chr>, inter2 <dbl>,
#   interaction <dbl>, civilian_targeting <chr>, iso <dbl>, region <chr>,
#   country <chr>, admin1 <chr>, admin2 <chr>, admin3 <chr>, location <chr>,
#   geo_precision <dbl>, source <chr>, source_scale <chr>, notes <chr>,
#   fatalities <dbl>, tags <chr>, timestamp <dbl>, population_1km <dbl>, …

First I create subsets for each event type, by default each event subset will inherit the SF object characteric

Battles <- filter(conflict_sf, event_type == "Battles")
Violence_CV <- filter(conflict_sf, event_type == "Violence against civilians")
Protests <- filter(conflict_sf, event_type == "Protests")
Riots <- filter(conflict_sf, event_type == "Riots")
Explosions <- filter(conflict_sf, event_type == "Explosions/Remote violence")
Strategic_dev <- filter(conflict_sf, event_type == "Strategic developments")
class(Battles)
[1] "sf"         "tbl_df"     "tbl"        "data.frame"

Visualising of Fatalities by Event Type

In the shiny App , we can enable users to filter by

  • specific year,

  • year range, and

  • event_type

  • range of number of fatalities

By using subsets of event types which have been converted to sf objects.

I will use the Geometry points from our data sets to plot the points of events in the maps from the leaflet package.

scaleFactor <- 2  

leaflet(Battles) %>%
  addProviderTiles(providers$CartoDB.Positron) %>%
  addPolygons(data = mmr_shp_mimu_1, color = "#444444", weight = 1, fillOpacity = 0.5) %>% # Adding borders
  
  addCircleMarkers(popup = ~paste("Event: Battles<br>State/Region:", admin1, 
                                  "<br>Actor1:", actor1, "<br>Actor2:", actor2,
                                  "<br>Year:", year, "<br>Fatalities:", fatalities),
                   radius = ~sqrt(fatalities) * scaleFactor,
                   fillColor = "red", fillOpacity = 0.4, color = "#FFFFFF", weight = 1) %>% 
  setView(lng = 96.1603, lat = 19.745, zoom = 6)
scaleFactor <- 2  

leaflet(Violence_CV) %>%
  addProviderTiles(providers$CartoDB.Positron) %>%
  addPolygons(data = mmr_shp_mimu_1, color = "#444444", weight = 1, fillOpacity = 0.5) %>% # Adding borders
  
  addCircleMarkers(popup = ~paste("Event: Violence on Civillians<br>State/Region:", admin1, 
                                  "<br>Actor1:", actor1, "<br>Actor2:", actor2,
                                  "<br>Year:", year, "<br>Fatalities:", fatalities),
                   radius = ~sqrt(fatalities) * scaleFactor,
                   fillColor = "red", fillOpacity = 0.4, color = "#FFFFFF", weight = 1) %>% 
  setView(lng = 96.1603, lat = 19.745, zoom = 6)
scaleFactor <- 2  

leaflet(Protests) %>%
  addProviderTiles(providers$CartoDB.Positron) %>%
  addPolygons(data = mmr_shp_mimu_1, color = "#444444", weight = 1, fillOpacity = 0.5) %>% # Adding borders
  
  addCircleMarkers(popup = ~paste("Event: Protests<br>State/Region:", admin1, 
                                  "<br>Actor1:", actor1, "<br>Actor2:", actor2,
                                  "<br>Year:", year, "<br>Fatalities:", fatalities),
                   radius = ~sqrt(fatalities) * scaleFactor,
fillColor = "red", fillOpacity = 0.4, color = "#FFFFFF", weight = 1) %>% 
  setView(lng = 96.1603, lat = 19.745, zoom = 6)
scaleFactor <- 2  

leaflet(Riots) %>%
  addProviderTiles(providers$CartoDB.Positron) %>%
  addPolygons(data = mmr_shp_mimu_1, color = "#444444", weight = 1, fillOpacity = 0.5) %>% # Adding borders
  
  addCircleMarkers(popup = ~paste("Event: Riots<br>State/Region:", admin1, 
                                  "<br>Actor1:", actor1, "<br>Actor2:", actor2,
                                  "<br>Year:", year, "<br>Fatalities:", fatalities),
                   radius = ~sqrt(fatalities) * scaleFactor,
                   fillColor = "red", fillOpacity = 0.4, color = "#FFFFFF", weight = 1) %>% 
  setView(lng = 96.1603, lat = 19.745, zoom = 6)
scaleFactor <- 2  

leaflet(Explosions) %>%
  addProviderTiles(providers$CartoDB.Positron) %>%
  addPolygons(data = mmr_shp_mimu_1, color = "#444444", weight = 1, fillOpacity = 0.5) %>% # Adding borders
  
  addCircleMarkers(popup = ~paste("Event: Explosions<br>State/Region:", admin1, 
                                  "<br>Actor1:", actor1, "<br>Actor2:", actor2,
                                  "<br>Year:", year, "<br>Fatalities:", fatalities),
                   radius = ~sqrt(fatalities) * scaleFactor,
fillColor = "red", fillOpacity = 0.4, color = "#FFFFFF", weight = 1) %>% 
  
  setView(lng = 96.1603, lat = 19.745, zoom = 6)